{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Example: Regenerating Data from\n",
    "# [R. Wu et al. / Elec Acta 54 25 (2010) 7394–7403](http://www.sciencedirect.com/science/article/pii/S0013468610009503)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy as sp\n",
    "import numpy as np\n",
    "import openpnm as op\n", "%config InlineBackend.figure_formats = ['svg']\n",
    "import matplotlib.pyplot as plt\n",
    "import openpnm.models.geometry as gm\n",
    "import openpnm.topotools as tt\n",
    "%matplotlib inline\n",
    "np.random.seed(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the workspace loglevel to not print anything"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws = op.Workspace()\n",
    "ws.settings[\"loglevel\"] = 50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run shared_funcs.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can run multiple times as the network sizes are randomly generated between a given range we can obtain an average"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_values = []\n",
    "y_values = []\n",
    "\n",
    "for ensemble in range(10):\n",
    "    x_ensemble, y_ensemble = simulation(n=8)\n",
    "    x_values.append(x_ensemble)\n",
    "    y_values.append(y_ensemble)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_values = np.asarray(x_values).flatten()\n",
    "y_values = np.asarray(y_values).flatten()\n",
    "plt.figure()\n",
    "from matplotlib.font_manager import FontProperties\n",
    "fontP = FontProperties()\n",
    "fontP.set_size('small')\n",
    "\n",
    "wu_average_x_values = [0.004, 0.021, 0.052, 0.081, 0.129, 0.162, 0.186, 0.219, 0.261,\n",
    "                       0.286, 0.324, 0.363, 0.42, 0.478, 0.531, 0.586, 0.64, 0.698, 0.747, 0.802]\n",
    "wu_average_y_values = [0.118, 0.113, 0.105, 0.096, 0.085, 0.078, 0.07, 0.062, 0.054, 0.049, 0.04,\n",
    "                       0.033, 0.027, 0.02, 0.012, 0.006, 0.003, 0.002, 0.002, 0.002]\n",
    "\n",
    "p1, = plt.plot(x_values, y_values, 'ko')\n",
    "p2, = plt.plot(wu_average_x_values, wu_average_y_values, 'ro')\n",
    "plt.title('normalized diffusivity versus saturation')\n",
    "plt.xlabel('saturation')\n",
    "plt.ylabel(r'$\\frac{D_e}{D_b}$')\n",
    "#plt.ylim([0, .15])\n",
    "plt.xlim([0, 1])\n",
    "plt.legend([p1, p2],\n",
    "                   [r'$\\frac{D_e}{D_b} = f(\\epsilon, \\phi)g(s, \\phi)$' + '\\n' + r'$X = 1.8$' +\n",
    "                   '\\n' + r'$Z_t = 2.0$' + '\\n' + r'$Z_i = 4.0$' + '\\n' + r'$\\beta = 1.0$' + '\\n' + r'$n = 14$', \"Wu's results\"])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally extract the g(S) function for relative diffusivity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "normalize_factor = max(y_values)\n",
    "g_values = y_values / normalize_factor\n",
    "\n",
    "wu_saturation = [0.004, 0.066, 0.0930, .119, 0.14, 0.175, 0.209, 0.24, 0.282, 0.32, 0.371, 0.413,\n",
    "                 0.464, 0.517, 0.605, 0.672, 0.761, 0.831, 0.898, 0.948, 0.996]\n",
    "wu_g_values = [0.986, 0.838, 0.758, 0.701, 0.651, 0.576, 0.516, 0.456, 0.39, 0.335, 0.268, 0.221,\n",
    "               0.171, 0.111, 0.067, 0.04, 0.019, 0.007, 0.003, 0.003, 0.003]\n",
    "\n",
    "p1, = plt.plot(x_values, g_values, 'ko')\n",
    "p2, = plt.plot(wu_saturation, wu_g_values, 'ro')\n",
    "plt.title('g(s) versus saturation')\n",
    "plt.xlabel('saturation')\n",
    "plt.ylabel('g(s)')\n",
    "plt.legend([p1, p2],\n",
    "                   [\"our values\", \"Wu's values (fitted curve)\"], loc='center left', bbox_to_anchor=(1, 0.5), prop = fontP)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
